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Abstract 


The project is based upon research of Tessler et al. (1994) on an improved variational 
formulation for post-processing stress predictions in Finite Element Analysis. The 
methodology, called Smoothing Element Analysis (SEA), employs a three-node 
smoothing finite element. The present effort focused on verifying the basic constant strain 
criterion for the three-node smoothing element subject to a set of internal penalty 
constraints. The convergence characteristics of the element are assessed by first deriving 
the constrained form of the assumed element stress and stress gradient fields, and then by 
verifying the validity of the constant strain criterion once the element penalty constraints 
are explicitly imposed. The analytical investigation is carried out with the use of the 
symbolic manipulation code Mathematica. 

Introduction 

An improved variational formulation called Smoothing Element Analysis (SEA), 
developed by Tessler et al., (1994), serves as a foundation for the enhancement of finite- 
element obtained deformation and stress response. In the case of stress predictions, C 1 - 
continuous stress field from a finite element solution is enforced into a C 1 -continuous 
stress field with continuous stress gradients. These enhanced results are ideally suited for 
error estimation since the stress gradients can be used to assess equilibrium satisfaction. 
The approach is employed as a post-processing step in finite element analysis. The 
variational statement combines the discrete-least squares, and penalty-constraint 
functionals, thus enabling automated recovery of smooth stresses and stress gradients. 

The practical issues whose adequate resolution is essential for a successful application 
of the approach are: 

(1) The SEA mesh and the number of the discrete stresses extracted from the Finite 
Element Analysis (FEA) mesh should be properly interrelated in order to produce a 
determined system of SEA equations. To fulfill this requirement, Tessler et al. (1994) 
proposed specific guidelines. An automated generation of the SEA mesh may also be 
necessary to make the post-processing transparent to the user. 

(2) The FEA stresses need to be extracted at the discrete elemental locations that are 
best suited for the recovery. Optimal (i.e., superconvergent) Barlow points and Gauss 
integration points have been successfully used. 

(3) The smoothing element should not exhibit locking — a pathological stiffening 
phenomenon commonly exhibited in penalty-constrained elements. In this connection, a 
judicious choice of the clement shape functions is key to avoiding locking. 

The purpose of this effort is investigate the influence of penalty constraints on the 
characteristics of the smoothing element used in Tessler et al. (1994). Particularly, the 
convergence characteristics of Tessler’s smoothing element are assessed by way of 
deriving the constrained form of the assumed element stress and stress gradient fields, and 
by verifying the validity of the constant strain criterion once the element penalty 
constraints are explicitly imposed. This analytical investigation is facilitated by the use of 
the symbolic manipulation code Mathematica. 
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Error Functional 


In this section an error functional proposed by Tessler et al. (1994) for a two- 
dimensional plane formulation is reviewed. It is assumed that within a two-dimensional 
region £2 = {xG^ 2 }, where x={jt,}, i= 1,2, represents a position vector in Cartesian 
coordinates, a finite element-derived stress field < 7 A (x) has been obtained by means of a 
discretization of £2 with characteristic element size h- The smoothed stress field, er'(x), 
is to be constructed from <r A (x) via a variational formulation. The variational statement 
involves scalar quantities only, and so each component of <J A (x) is smoothed 
independently. Hence, in the following reference is made only to components <7* and c f . 
The finite element stress field is sampled at x q , q = 1,2, ..., N, to obtain the set of 
stresses {o^}, i.e., =cr*(x ? ). The sampled stresses are those extracted at the Gauss 

integration points, Barlow points, or other element locations in the finite element 
analysis. To minimize the error functional, we adopt the finite element methodology and 
therefore discretize £2 with n e/ “recovery” or “smoothing” finite elements such that 
£2 = u"*,£2% where £2 f is the domain of smoothing element e Within our recovery element 
model, we use C°-continuous interpolation functions for the stress, <f , and the 
independent quantities 9- , 7=1, 2, whose mathematical interpretation will be readily 
established. The error functional to be minimized can be written as 

<i>= ^X w K- <t, K)] 2+ 

where w q and p(x) are the appropriate weight functions; y is a normalization factor; X is 

a dimensionless parameter; and a comma denotes partial differentiation. Because the 
highest partial derivative in equation ( 1 ) is of order one, the field variables need only be 
approximated with C°-continuous shape functions. 

The first term in equation (1) represents a discrete least-squares functional in which 
the squared ‘error’ between the smoothed stress field and the sample data is computed for 
all sampled stresses. The term can be normalized in several different ways; presently, the 
normalization factor equals the total number of the sampled stresses, i.e., y=N . The 
discrete weights w q are introduced so that sample data known to be of higher accuracy 

can be assigned more weight than less accurate data. 

The second term in equation (1) represents a penalty functional which, for X 
sufficiently large, enforces the derivatives of the smoothed stress field & . to approach the 

corresponding 9* variable pointwise, i.e., 

(z'-x,y) in £2 f (2) 


Jp(x)[k -a;) 2 +k -0;> 2 ]^ (1) 
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Theoretically, the greater the value of A, the closer the correlation between <f and 9-, 

where C 1 continuity of & is achieved as A -> °o. In practice, A needs to be sufficiently 
large in order to enforce conditions (2), yet it should not be excessively large to cause ill- 
conditioning of the smooth solution. Because 6f are interpolated with continuous 
functions, the smoothed stress field, for all practical purposes, is C 1 continuous. The 
weight function p(x) is introduced in the functional to allow the enforcement of C 1 
continuity to be somewhat relaxed in certain regions of ft and more strictly enforced in 
others, if so desired. Also note that by specifying the weights w q and p(x) to vanish in 

regions outside a given domain of interest, a local or patch analysis is admitted. Presently, 
we only consider the special case where w q = 1 and p(x) = 1 , that is all stress data are 
treated equally and the C 1 continuity is enforced throughout the Q domain. 

Assumed Element Fields 

Although the functional (1) admits C°-continuous shape functions for the field 
variables (f, (%, and 0 s ,, the constraints (2) impose certain restrictions on the suitable 
choice of shape functions. In a similar plate theory formulation, constraints of this type 
are known to cause locking (i.e., severe stiffening) when conventional isoparametric 
interpolations are used. (In the present context, locking would manifest itself in a 
smoothed stress field, <f , that grossly underestimates the ‘true’ stress distribution.) 
When Cf is interpolated with a polynomial one degree higher than those for the 9\ and 
% variables, using anisoparametric interpolations, the locking effect is alleviated or 
completely eliminated (Tessler, 1985). 

Another important consideration is the nodal configuration that is best suited for the 
smoothing element. It turns out that a three-node triangle is well-suited for this purpose 
because (a) from a modeling standpoint, it represents the most versatile element topology, 
and (b) it permits a one-to-one linear mapping between the global and element local (area- 
parametric) coordinates, thus allowing a straightforward identification of the sampled 
stress data within the smoothing element. 

The anisoparametric interpolations for a three-node element involve quadratic 
approximation of <f and linear approximations of ff x and 6? which can be expressed in 
matrix form as 

o s ~ za e + m 9' + 10;, 0; = z 0' (i=l,2) (3) 

where <f , ff. are 3x1 vectors of nodal degrees-of-ffeedom (dof), z is selected as a row- 
vector of a linear shape function, and m and 1 are selected as row-vectors of quadratic 
shape functions. Their explicit forms given in terms of area-parametric coordinates are 

z = {z,,z 2 ,z 3 }, m = {m,,m 2 ,m 3 }, ! = {/,, / 2 ,/ 3 ) (3.1) 

where 
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Zi = 2 t(c, + V + a,y), m, = 4(<W; -^z.-z*), /, = t(*W* -&*z,-z>) 
a , = ** -■*>» b t = y } - y *> <■*,• = jc y y t -^ t y y - (/ = 1 , 2 , 3 , j = 2 , 3 , 1 , £ = 3 , 1 , 2 ) 

and A denotes the area of a triangular element. 

Note that these interpolations are consistent with a three-node element which has only 
three dof per node, even though <f is quadratic and 0* are linear functions. Moreover, 
equations (3) ensure that the gradient of the smoothed stress, <f n is the same degree 
polynomial as that representing 6f, i.e., they are both linearly distributed across . This 
naturally leads to a reasonable expectation that penalty constraints (2) can be adequately 
fulfilled without over constraining ( locking ) the element. 

Edge Penalty Constraints 

A straightforward manipulation of the two constraint equations in (2), in which 
equations (3) are introduced, produces three edge-wise constraints per element. For the 
element edge defined by nodes i and j, the edge constraint equation has a simple form in 
terms of the nodal dof corresponding to the edge (Pomeranz, 1995; Tessler, 1985): 


a; -a; -> j(x ( . -Xjxer +e;,)+{(y l - yj w yi , + 0 ;.) 


(4) 


where x k and y k (k = ij ) are the nodal coordinates. 

The three edge constraint equations ensure that there are only six independent dof per 
element, thus properly describing the complete parabolic field of O'. They also facilitate a 
simple calculation of the total number of independent dof in the mesh. The key aspect of 
these constraints is that they control the mechanisms of locking. Their assessment in the 
context of assembly of elements can provide proper insight into preferable discretization 
patterns for such elements. For example, a fully non-locking behavior is achieved by 
producing SEA meshes made of quadrilateral macro-elements that are formed with four 
triangles in a cross-diagonal pattern. 

Constant Strain Criterion 


Let us consider an arbitrary triangular element as shown in the diagram below. 
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The satisfaction of the constant strain criterion in the finite element method ensures 
convergence of the method as the mesh is refined. In this case, it is expected that each 
individual finite element accommodates constant strains. Mathematically, the criterion is 
verified by summing up on all shape functions for each field that is approximated, and the 
resulting sum should add up to unity. This can be readily verified for the unconstrained 
element fields in (3). 

The constraint equations for the edges of this element can be written as follows. 


Edge 1-2 

- a 1 = i ((*, - *2 )( 0 X , + e x 2 ) + (y, - y 2 )(0 vl + e y2 )) 

Edge 2-3 

°2 = f ((* 2 -* 3 X 0*2 +0* 3 ) + (t 2 + 0 , 3 )) 

Edge 3-1 

- G \ - 2 ((*3 -*1X0*3 + 0*l) + (T3 X0,3 + 0 y , )) 

(5) 

Using Mathematica, the three constraint equations are solved for <r,, a 2 , and 6 xi . 
When these solutions are substituted into the original definitions for a s , G 5 X , and 0 ¥ \ the 
following expressions of the three element fields are derived 

=^1^3 +£ 2 (*)0*i +g 3 (*)0* 2 +g 4 (**>')0vi +g 5 (*,y)0 y2 +8 6 (x,y)e y3 

0* =d t (x)6 xl + d 2 (x)O x2 + d 2 (x,y)6 yl + d 4 (x,y)0 y2 +d 5 (x,y)O yi ( 6 ) 

0y =e x (*,y)0 vl + e l (x,y)0 y2 +e i (x,y)0 vi 


where g, d, and e are shape functions whose expressions are summarized in the Appendix. 

To verify the constant strain criterion for the resulting element fields, the summation of 
the g, d, and e shape functions is carried out with the use of Mathematica. The resulting 
equations are as follows 

S = = l + x-x i +y~y i , </ = 2X=l, e = ^e,=l (7) 

'=1 t=l i=l 

Note that both 0 x and 9 y fulfill the constant strain criterion fora finite size element since d 

=1 and e=\. On the other hand, g only approaches unity in the limit as the element size 
diminishes to zero, i.e., 

* -» * 3 > y -> y 3 ( 8 ) 
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giving rise to g -» 1 . Thus, the element convergence is ensured as the smoothing mesh is 
refined. 

Conclusions 

The project has been a continuation of the research of Tessler ct al. (1994) on an 
improved variational formulation for post-processing stress predictions in Finite Element 
Analysis. The effort focused on verifying the basic constant strain criterion for the three- 
node smoothing element subject to a set of internal penalty constraints. The convergence 
characteristics of the element were assessed by first deriving the constrained form of the 
assumed element stress and stress gradient fields through a process of simplifications 
using Mathematica. The element penalty constraints were explicitly imposed using the 
formulas set out by Tessler et al. (1994). Then the validity of the constant strain criterion 
for the element was verified. As the smoothing mesh is refined, the constant strain 
criterion is satisfied. 
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Appendix 


The following are the shape functions in (6) as solved by Mathematica™. 
g i= 1 

g 2 = ((-x + x3)*(x - 2*x2 + x3))/(2*(-xl + x2)) 
g 3 = ((x - x3)*(-x + 2*xl - x3))/(2*(xl - x2)) 


g 4 = (xl*x2*y A 2 - x2 A 2*y A 2 - xl*x3*y A 2 + x2*x3*y A 2 - 2*x*xl*y*y2 + 2*x*x2*y*y2 + 
2*xl*x3*y*y2 - 2*x2*x3*y*y2 + x A 2*yl*y2 - 2*x*x2*yl*y2 + 2*x2*x3*yl*y2 - 
x3 A 2*yl*y2 - x A 2*y2 A 2 + 2*x*xl*y2 A 2 - 2*xl*x3*y2 A 2 + x3 A 2*y2 A 2 + 2*x*xl*y*y3 - 
2*x*x2*y*y3 - 2*xl*x2*y*y3 + 2*x2 A 2*y*y3 - x A 2*yl*y3 + 2*x*x2*yl*y3 - 
2*x2*x3*yl*y3 + x3 A 2*yl*y3 + x A 2*y2*y3 - 2*x*xl*y2*y3 + 2*xl*x3*y2*y3 - 
x3 A 2*y2*y3 + xl*x2*y3 A 2 - x2 A 2*y3 A 2 - xl*x3*y3 A 2 + x2*x3*y3 A 2)/(2*(-xl + x2)*(- 
(x2*yl) + x3*yl + xl*y2 -x3*y2 -xl*y3 + x2*y3)) 

gs = (-(xl A 2*y A 2) + xl*x2*y A 2 + xl*x3*y A 2 - x2*x3*y A 2 + 2*x*xl*y*yl - 2*x*x2*y*yl 
- 2*xl*x3*y*yl + 2*x2*x3*y*yl - x A 2*yl A 2 + 2*x*x2*yl A 2 - 2*x2*x3*yl A 2 + 
x3 A 2*yl A 2 + x A 2*yl*y2 - 2*x*xl*yl*y2 + 2*xl*x3*yl*y2 - x3 A 2*yl*y2 - 
2*x*xl*y*y3 + 2*xl A 2*y*y3 + 2*x*x2*y*y3 - 2*xl*x2*y*y3 + x A 2*yl*y3 - 
2*x*x2*yl*y3 + 2*x2*x3*y]*y3 - x3 A 2*yl*y3 - x A 2*y2*y3 + 2*x*xl*y2*y3 - 
2*xl*x3*y2*y3 + x3 A 2*y2*y3 - xl A 2*y3 A 2 + xl*x2*y3 A 2 + xl*x3*y3 A 2 - 
x2*x3*y3 A 2)/(2*(-xl + x2)*(-(x2*yl) + x3*yl + xl*y2 - x3*y2 - xl*y3 + x2*y3)) 

g6= (-(xl*y) + x3*y + x*yl - x3*yl - x*y3 + xl*y3)/(2*(-xl + x2)) + (x2*y - x3*y - 
x*y2 + x3*y2 + x*y3 - x2*y3)/(2*(-xl + x2)) + ((-x + x3)*(yl - y2)*(-(xl*y) + x2*y + 
x*yl - x2*yl - x*y2 + xl*y2))/(2*(-xl + x2)*(x2*yl -x3*yl -xl*y2 + x3*y2 + xl*y3 - 
x2*y3)) + ((xl*y - x2*y - x*yl + x2*yl + x*y2 - xl*y2)*(y - y3))/(2*(x2*yl - x3*yl - 
xl*y2 + x3*y2 + xl*y3 -x2*y3)) 

d x = (x - x2)/(xl - x2) 

d 2 = (x - xl)/(-xl + x2) 

^ 3 = ((-(xl*y) + x2*y + x*yl -x2*yl - x*y2 + xl*y2)*(y2 - y3))/((-xl + x2)*(-(x2*yl) + 
x3*yl + xl*y2 - x3*y2 - xl*y3 + x2*y3)) 

d 4 = ((xl*y - x2*y - x*yl + x2*yl + x*y2 - xl*y2)*(-yl + y3))/((-xl + x2)*(x2*yl - 
x3*yl - xl *y2 + x3*y2 + xl *y3 - x2*y3)) 
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ds= ((y I - y2)*(-(xl*y) + x2*y + x*yl -x2*yl - x*y2 + xl*y2))/((-xl + x2)*(-(x2*yl) + 
x3*yl + xl *y2 - x3*y2 - xl *y3 + x2*y3)) 

e\= (x2*y - x3*y - x*y2 + x3*y2 + x*y3 - x2*y3)/(x2*yl - x3*yl - xl*y2 + x3*y2 + 
xl*y3 - x2*y3) 

e 2 = (xl*y - x3*y - x*yl + x3*yl + x*y3 - xl *y3)/(-(x2*yl) + x 3* y i + x i* y 2 . x 3* y 2 - 
xl*y3 + x2*y3) 

ey= (-(xl *y) + x2*y + x*yl -x2*yl - x*y2 + xl *y2)/(-(x2*yl) + x3*yl + xl*y2 - x3*y2 
- xl*y3 + x2*y3) 


Notation 

In the above expressions, xl=jc,, yl=^,, xl A 2=jc t 2 , and the asterisk (*) denotes 
multiplication. 
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